using EMTSim
using BlockSystems
using ModelingToolkit
using NetworkDynamics
using Graphs
using OrdinaryDiffEq
using DiffEqCallbacks
using SteadyStateDiffEq
using Plots
using Unitful
using CSV
using DataFrames

Constants and unit stuff.

ω0    = 2π*50u"rad/s"
Sbase = 300u"MW"
Vbase = 110u"kV" #* sqrt(2/3)
Ibase = Sbase/(Vbase )#* √3) # why the √3 ?
Cbase = Ibase/Vbase
Lbase = Vbase/Ibase
Rbase = (Vbase^2)/Sbase

Rline = 1u"Ω" / Rbase           |> u"pu"
Pload = 300u"MW" / Sbase        |> u"pu"
Rload = (1u"pu")^2 / Pload      |> u"pu" # R=U^2/P
Cline = (2e-6)u"F" / Cbase      |> u"s"
Lline = (1/100π)u"H" / Lbase    |> u"s"

Slack Bus

@variables t u_d(t) u_q(t)
dt = Differential(t)

slackblock = IOBlock([dt(u_d) ~ 0, dt(u_q) ~ 0], [], [u_d, u_q]; name=:slack)

IOBlock :slack with 2 eqs
  ├ inputs:  (empty)
  ├ outputs: u_d(t), u_q(t)
  ├ istates: (empty)
  └ iparams: (empty)

create ODE Vertex from this block

slack = ODEVertex(slackblock)

R Load

@variables t i_d(t) i_q(t)
@parameters u_d(t) u_q(t) R
loadblock = IOBlock([i_d ~ -1/R * u_d,
                     i_q ~ -1/R * u_q],
                    [u_d, u_q], [i_d, i_q],
                    name=:load)

IOBlock :load with 2 eqs
  ├ inputs:  u_d(t), u_q(t)
  ├ outputs: i_d(t), i_q(t)
  ├ istates: (empty)
  └ iparams: R

The load is used as a current source in a BusBar

busblock = BusBar(loadblock; name=:loadbus)
busblock = set_p(busblock, Dict(:C=>ustrip(u"s", Cline), :ω0=>ustrip(u"rad/s", ω0)))

IOBlock :loadbus with 2 eqs
  ├ inputs:  i_d(t), i_q(t)
  ├ outputs: u_d(t), u_q(t)
  ├ istates: (empty)
  └ iparams: load₊R

The BusBar block can be used to create an ODE Vertex

load = ODEVertex(busblock, [:load₊R])

ODE edge for the RL Line

@variables t i_d(t) i_q(t)
@parameters u_d_src(t) u_q_src(t) u_d_dst(t) u_q_dst(t) R L ω

lineblock = IOBlock([dt(i_d) ~  ω * i_q  - R/L * i_d + 1/L*(u_d_src - u_d_dst),
                     dt(i_q) ~ -ω * i_d  - R/L * i_q + 1/L*(u_q_src - u_q_dst)],
                    [u_d_src, u_q_src, u_d_dst, u_q_dst],
                    [i_d, i_q],
                    name=:RLLine)
lineblock = set_p(lineblock, Dict(:R=>NoUnits(Rline), :L=>ustrip(u"s", Lline), :ω=>ustrip(u"rad/s", ω0)))

IOBlock :RLLine with 2 eqs
  ├ inputs:  u_d_src(t), u_q_src(t), u_d_dst(t), u_q_dst(t)
  ├ outputs: i_d(t), i_q(t)
  ├ istates: (empty)
  └ iparams: (empty)

This block can be used to create an ODEEdge

edge = ODEEdge(lineblock)

Set up the network

g = SimpleGraph(2)
add_edge!(g, 1, 2)
nd = network_dynamics([slack, load], edge, g)

we use SteadyStateDiffeq to find the steady state

uguess = zeros(length(nd.syms))
uguess[1] = 1.0 # set the d component of the slack to 1 from the beginning
p = ([0, NoUnits(Rload)], nothing)
ssprob = SteadyStateProblem(nd, uguess, p)
u0 = solve(ssprob, DynamicSS(AutoTsit5(Rosenbrock23())))

# Simulation

We simulate a disconnection of the load at t=0.1 s

tspan = (0.0, 0.2)
cb = PresetTimeCallback(0.1, integrator -> integrator.p = ([0, Inf], nothing))
prob = ODEProblem(nd, u0, tspan, p; callback=cb)
sol = solve(prob, Tsit5())

# plot the dq voltage of node 2
# plot(sol, vars=[3,4]; xlims=(0.095, 0.124))

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 1331-element Vector{Float64}:
 0.0
 0.004596192622990568
 0.0048673680168831515
 0.005112978790178177
 0.005440094588913402
 0.005671217534552417
 0.0059000150078427115
 0.006246626451068669
 0.006552133024265822
 0.006826198260014092
 ⋮
 0.1989440189891149
 0.19909314942295595
 0.19923127797379012
 0.1993984607587681
 0.1995310996834446
 0.1996784230803493
 0.19983249459047903
 0.19996739423890775
 0.2
u: 1331-element Vector{Vector{Float64}}:
 [1.0, -2.8853683219220016e-22, 0.9758038843594189, -0.024221291204110477, 0.9764177046988668, 0.000507707199198463, -0.9764177046988668, -0.000507707199198463]
 [1.0, -2.8853683219220016e-22, 0.9757986924082586, -0.02422368266229332, 0.9764130606766801, 0.0005111425262533142, -0.9764130606766801, -0.0005111425262533142]
 [1.0, -2.8853683219220016e-22, 0.9758024273244535, -0.024226587122023014, 0.9764123064894793, 0.0005064159745894732, -0.9764123064894793, -0.0005064159745894732]
 [1.0, -2.8853683219220016e-22, 0.9758052654054002, -0.02422284293245541, 0.9764170035367953, 0.0005057270922529847, -0.9764170035367953, -0.0005057270922529847]
 [1.0, -2.8853683219220016e-22, 0.9758117316972282, -0.024222076573326895, 0.976420833145115, 0.0005003714301078855, -0.976420833145115, -0.0005003714301078855]
 [1.0, -2.8853683219220016e-22, 0.9758045116236415, -0.024220187042279163, 0.9764189884120059, 0.0005076895743561899, -0.9764189884120059, -0.0005076895743561899]
 [1.0, -2.8853683219220016e-22, 0.9758037781083057, -0.02422116133412916, 0.9764177679573095, 0.0005078646034726465, -0.9764177679573095, -0.0005078646034726465]
 [1.0, -2.8853683219220016e-22, 0.9758029539373411, -0.024221133678456772, 0.9764173908031514, 0.0005086081298502008, -0.9764173908031514, -0.0005086081298502008]
 [1.0, -2.8853683219220016e-22, 0.9758020633459226, -0.024222572316079737, 0.9764156824341914, 0.0005086959362062117, -0.9764156824341914, -0.0005086959362062117]
 [1.0, -2.8853683219220016e-22, 0.9758036600018949, -0.0242236111649951, 0.976415540326711, 0.0005067753622784793, -0.976415540326711, -0.0005067753622784793]
 ⋮
 [1.0, -2.8853683219220016e-22, 1.000669569462355, -0.0006422832059868946, -1.2912966491950219e-5, 0.02531842394222369, 1.2912966491950219e-5, -0.02531842394222369]
 [1.0, -2.8853683219220016e-22, 1.0005880351616798, -0.0006608789007876454, -1.3711426617890149e-5, 0.025384378994485632, 1.3711426617890149e-5, -0.025384378994485632]
 [1.0, -2.8853683219220016e-22, 1.0006069611413242, -0.0005980490359880196, 6.192248318034841e-5, 0.02538295439472409, -6.192248318034841e-5, -0.02538295439472409]
 [1.0, -2.8853683219220016e-22, 1.000679130244612, -0.0006254352872959244, 9.226484959078049e-6, 0.025317823783546686, -9.226484959078049e-6, -0.025317823783546686]
 [1.0, -2.8853683219220016e-22, 1.0006160890452325, -0.0006681106424376321, -3.3780972436427184e-5, 0.02536076536537274, 3.3780972436427184e-5, -0.02536076536537274]
 [1.0, -2.8853683219220016e-22, 1.0005851144920301, -0.0006143909445733833, 4.3558069089972444e-5, 0.02539384688969939, -4.3558069089972444e-5, -0.02539384688969939]
 [1.0, -2.8853683219220016e-22, 1.0006698769905904, -0.00060301246168973, 4.5377922772656195e-5, 0.025329643442160883, -4.5377922772656195e-5, -0.025329643442160883]
 [1.0, -2.8853683219220016e-22, 1.0006510657449301, -0.0006605567589685328, -3.0054199486610306e-5, 0.02533793333598354, 3.0054199486610306e-5, -0.02533793333598354]
 [1.0, -2.8853683219220016e-22, 1.000630835548742, -0.0006658911685662713, -3.497205532094123e-5, 0.025352879465425606, 3.497205532094123e-5, -0.025352879465425606]

Transformation of the results back to a,b,c frame

a,b,c = Tdqinv(sol.t, sol[3,:], sol[4,:])

Comparison with power factory results

Read the Power Facory data and plot for reference

$V_{base}$ is RMS phase-phase voltage.

\[V_{star} = \frac{1}{\sqrt{3}} * V_{base}\\ \hat{V} = \sqrt{\frac{2}{3}} * V_{base}\]

Therefore, i have to multiply the a, b and c results by sqrt(3/2)

df = CSV.read(joinpath(dirname(pathof(EMTSim)), "..", "data","PowerFactory", "Test_EMT.csv"), DataFrame, skipto=3,
              header=[:t, :u_1_a, :u_1_b, :u_1_c, :u_2_a, :u_2_b, :u_2_c])

xlims = (0.099,0.124)
plot(sol.t, a.*sqrt(3/2); label="u_2_a", xlims)
plot!(df.t, df.u_2_a; label="PowerFactory A")

plot!(sol.t, b.*sqrt(3/2); label="u_2_b")
plot!(df.t, df.u_2_b; label="PowerFactory B")

plot!(sol.t, c.*sqrt(3/2); label="u_2_c")
plot!(df.t, df.u_2_c; label="PowerFactory C")

... lets have a closer look

xlims!(0.0995,0.105)

This page was generated using Literate.jl.